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ABSTRACT 

We create dynamical models of the massive elliptical galaxy, NGC 4649, using the 
N-body made-to-measure code, NMAGIC, and kinematic constraints from long-slit 
and planetary nebula (PN) data. We explore a range of potentials based on previous 
determinations from X-ray observations and a dynamical model fitting globular cluster 
(GC) velocities and a stellar density profile. The X-ray mass distributions are similar 
in the central region but have varying outer slopes, while the GC mass profile is higher 
in the central region and on the upper end of the range further out. Our models cannot 
differentiate between the potentials in the central region, and therefore if non-thermal 
pressures or multi-phase components are present in the hot gas, they must be smaller 
than previously inferred. In the halo, we find that the PN velocities are sensitive tracers 
of the mass, preferring a less massive halo than that derived from the GC mass profile, 
but similar to one of the mass distributions derived from X-rays. Our results show that 
the GCs may form a dynamically distinct system, and that the properties of the hot 
gas derived from X-rays in the outer halo have considerable uncertainties that need 
to be better understood. Estimating the mass in stars using photometric information 
and a stellar population mass-to-light ratio, we infer a dark matter mass fraction in 
NGC 4649 of - 0.39 at IR^ (10.5 kpc) and - 0.78 at AR^. We find that the stellar 
orbits are isotropic to mildly radial in the central ^ 6 kpc depending on the potential 
assumed. Further out, the orbital structure becomes slightly more radial along R and 
more isotropic along z, regardless of the potential assumed. In the equatorial plane, 
azimuthal velocity dispersions dominate over meridional velocity dispersions, implying 
that meridional velocity anisotropy is the mechanism for flattening the stellar system. 

Key words: galaxies: ellipticals and lenticulars, CD - galaxies: kinematics and dy- 
namics - galaxies: evolution - galaxies: individual: NGC 4649 (M60) - planetary neb- 
ulae: general - dark matter 



1 INTRODUCTION 

Massive elliptical galaxies are huge conglomerates of stars, 
dark matter and hot gas, often residing at the centres of 
dense environments. They are believed to evolve through a 
complex formation process, manifested in the intricate struc- 
ture of their orbits. Understanding the distribution of mass 
in massive ellipticals is vital for obtaining constraints on the 
dark matter content and on the orbital structure of the stars. 

The hot, low-density gas surrounding massive elliptical 
galaxies produces X-ray spectra consisting of emission lines 
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and continuous emission from thermal bremsstrahlung ra- 
diation. The spectra can be modelled to obtain the density 
and temperature profiles of the gas. If the gas is relatively 
undisturbed then we can assume hydrostatic equilibrium 
and derive the total mass distribu tion from the den- 
sity an d temperature p rofiles ( e.g. Nulsen &: Bohringeij 



1 19951; iFukazawa et al.1 l2006l: 'Humphrey et al.' '20061: 



i Churazov et al.„2008i : iNagino &: Matsushita,2009. : ,Das et al] 

boid i. 

Another common method for obtaining mass pro- 
files of massive elliptical galaxies is through the cre- 
ation of dynamical models. They can be const r ucted 
by superposing a library of orbits (e.g. iRix et al] Il997l : 



p. Das et al. 



Gebhardt et al.ll2003l: iThomas et ai]|2004l : [Cappellari et al] 
20061 : Ivan den Bosch et al.l |2008|) or distribution func- 



tions (e.g. Deiongh eet al ] Il996l : iGerhard et al ] 1 19981 : 
iKronawitter et al.. ,20Q(]| ), or by constructing a syst em of 
particles (e.g. NMAGIC, Ide Lorenzi et al.ll2008l . |2009| ) such 
that the projection of the system best reproduces observed 
surface-brightness and kinematic profiles. These give the 
mass distribution and orbital structure of the galaxies si- 
multaneously. 

Obtaining the mass distribution from X-rays is obvi- 
ously a simpler way and this can then be used as input in 
the dynamical models of galaxies, therefore mitigating the 
usual mass-anisotropy degeneracy. This would provide more 
stringent constraints on the orbital structure derived from 
dynamical models. Before doing this however, one has to 
understand why there are discrepancies between mass pro- 
files determined from X-rays and from d ynamical models , 
and how significant these discrepancies are. lDas et al.l (|2010l ) 
compared their X-ray mass determinations for a sample of 
six nearby massive elliptical galaxies to the most recent dy- 
namical mass determinations with inconclusive results. The 
most common discrepancy between the two is in the central 
~ 10 kpc. The X-ray circular velocity curve, 
is on average about 20% lower th an the dynamica l resul t 
in this re g ion. B ased on work by IChurazov et all HoiO), 
iDas et all (|2010l ) speculated that the most likely sources of 
this discrepancy are: 1) Non-thermal contributions to the 
pressure that are not considered in the application of hydro- 
static equilibrium, 2) multi-phase components in the gas, 3) 
a lack of spatially extended observational constraints in the 
dynamical models, and 4) mass profiles in the dynamical 
models that are not sufficiently general. 

In this paper, we would like to look in more detail 
at NGC 4 649, one of the m assive elliptical galaxies in the 
sample of iDas et all (|2010l ). residing at the centre of a 
sub-clump in the Virgo cluster. It is the fourth brightest 
early-type galaxy in th e clus ter. Long-slit k i nemat ics from 
iDe Bruvne et al] l|200lh and IPinknev et"aLl l|2003f ) show a 
high velocity dispersion of ~ 400 km/s in the centre, and a 
mean velocity of ~ 100 km/s along the major axis. There is 
also the rece nt catalogue of 121 GC line-of-sight (LOS) ve- 
lociti es fromlHwang et al] (|2008l ). building on the catalogue 
from [Bridges et aL 1 20061 ). They found significant rotation 
in the GC system, of order 141 km/s, and an average velocity 
dispersion of 234 km/s. 

There are several mass distributions deter- 
mined for N GC 4649 from X- ray data using 
ROSAT data dTrinchieri et al.l 119971). Chandr a data 
llBrighenti fc Mathewd Il997l : iHumphrev et al.l |2006| . 
l2008l ). and a combination of Chandra and XMM-Newton 
data llNagino fc Matsushital l2009l : IChurazov et al.l l2010l : 
iDas et aLlboiOl) . The mass profiles ar e simila r in th e central 
~ 12 kpc except that of JHumphrex et al.1 (|2006i ). which 
is a little higher. Further out they all point towards the 
existence of a massive dark matter halo, but the precise 
value of th e circular velocity varies between ^ 380-500 km/s 
at 20 kpc (|Hwang et al.ll2008l : lOas et al.ll201Gl ). 

There are also several dynamical models in the lit- 
erature combinin g the photometric and kinematic data 
described above. [Bridges et al. |"2006l) u sed a xisymmetric 
Schwarzschild models and Hwang et al.l l|2008l ) used Jeans 
models to fit GC LOS velocities and long-slit kinematics 



in the X-ray potential from iHumphrev et al.1 (|2006l ). They 
found an isotropic t o a m odestly tangential orbital structure. 
IShen fc Gebhardd (|201G| ) used axisymmetric Schwarzschild 
models to fit additionally kinematics from the Hubble Space 
Telescope (HST) and carried out an independent mass anal- 
ysis. They confirmed the presence of a supermassive black 
hole in the centre and a dark matter halo. Their best-fit 
mass profile is higher than the mass profiles derived from 
X-ray observations in the central ~ 12 kpc, and i n the outer 
parts agrees best with the mass distribution of iDas et al.l 
|20l3). 

In this work, we use the made-to-mea sure particle- 
based code NMAGIC (|de Lorenzi et al.|[2007l ) to create dy- 
nami cal models of NGC 464 9. We use photometric profiles 
from iKormendv et al.l (|2009D . long-slit kinematic data from 
iPinknev et al.1 (|2003| ). and LOS velocities measured for a 
new catalogue of 298 PNe described in a companion pa- 
per (|Teodorescu et al.ll201ll ). The PNe trace the kinematics 
out to about 416", similar to the GCs, but with more than 
double the number of tracers. There is also strong evidence 
that PNe are good tracers of the stellar density and kinemat- 
ics in early-type galaxies (jCoccato et al.l[2009l ). We explore 
mass d istri butions based on the det erminations of lDas et al.l 
(|2010l ) and lShen fc Gebhardtl (|201G| ). which only differ in the 
central ~ 12 kpc. 

With our models we wish to address the following ques- 
tions: 

(i) How massive is the dark matter halo in NGC 4649? 

(ii) What is the orbital structure of the stars in NGC 
4649? 

(iii) Is the potential derived from X-ray observations ac- 
curate enough to determine dark matter mass fractions and 
the orbital structure in the halo of massive elliptical galax- 
ies? 

(iv) Are the dark matter mass fractions and orbital struc- 
ture in the halo derived from PNe consistent with that de- 
rived from GCs? 

We assume that N GC 4649 is at a distance of 16.83 Mpc 
(|Tonrv et al.1 12001^ and therefore 1" = 82 pc and 1 kpc 
= 12". We also assum e an effective radius, = 10.5 kpc 
(|Kormendv et aLll2009l ). Intrinsic radii are denoted by r and 
quoted in kpc, and projected radii are denoted by R and 
quoted in arcsec. Radii in cylindrical coordinates are also 
denoted by R. 

In Sections [5] and |31 we describe the photometric and 
kinematics constraints we use for the projected distribution 
function of the stars. In Section [4] we describe the dynami- 
cal models we create with NMAGIC and discuss the implica- 
tions of our results in Section[5l We end with our conclusions 
in Section [6l 



2 CONSTRAINTS ON THE DISTRIBUTION 
OF STARS 

2.1 Photometric data 



We use the l/-band photometry of iKormendv et"ai] (|2009| ) 
that combines new measurements with published profiles, 
and extends out to 693" along the major axis. Figure [T] 
shows the measured surface-brightness and ellipticity pro- 
files. The surface-brightness profile has a central core with a 
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Figure 1. Projected light distribution in NGC 4649: The black 
circles in the top panel show the measured V-band surface- 
brightness profile and the lines show the reproj acted surface- 
brightness profiles from the axisymmetric deprojections. The 
black circles in the bottom panel shows the measured elliptic- 
ity profile and the lines show the ellipticity profiles of the repro- 
jected light distributions from the axisymmetric deprojections. 
We carry out deprojections assuming inclinations of 45° (solid), 
60° (dashed), 75° (dotted), and 90° (dash-dotted). 



shallow decay ext ending to about 4" and then falls off more 
steeply outwards. iKormendv et al.1 (I2OO9I ) fit a Sersic profile 
between 5-488", finding a Sersic index, n = 5.36. The ellip- 
ticity of the isophotes is around 0.1 in the central 2.5" but 
then increases to about 0.2 in the outer parts. The average 
position angle (measured from North tow ards East in the 
sky) o f the major axis from the profile in iKormendv et al] 
hOO^ is 102° ±6° , consistent with the value of 105° adopted 
bv lPinknev et all 42003) for the kinematic data (see Section 
13. ip . Therefore we also adopt a value of 105° for the major 
axis. 
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Figure 2. Kinematic constraint s for NGC 4649: The fill ed circles 
show the positions of the PNe l|Teodorescu et al.|[2oTlh and the 
solid lines show the positions of the slits. The x-axis is along the 
major axis of the galaxy and the j/-axis is along the minor axis. 
Directions of North and East are shown in the top-left corner. 
One effective radius (-Re) is illustrated by the dashed, black line. 



inclination, 



2 2 . , ,2 . 2 . 

q — cos t + t, sm i 



(1) 



The maximum apparent flattening from the isophotal anal- 
ysis is 0.79 and if we assume a maximum intrinsic flattening 
of 0.5 then the most face-on inclination allowed is 45°. 

Figure [T] shows the reprojection of the output depro- 
jected profiles compared with the original input surface- 
brightness and ellipticity profiles. They fit the measured pro- 
files very well except in the central arcsec. Here the ellipticity 
profile is not so well defined, because the isophotes are close 
to circular as a result of seeing. 



2.2 3-D density distribution of the stars 

To obtain the intrinsic distribution of the stars we assume 
that th ey form an oblate , axisymmetric system. We use the 
code of lMagorrianI (Il999l ) that finds a smooth axisymmetric 
density distribution consistent with the surface-brightness 
distribution, for some assumed inclination angle. It chooses 
the solution that maximises a penalised likelihood, ensuring 
a smooth 3-D luminosity density. 

We carry out axisymmetric deprojections for inclina- 
tions of i = 45° , i = 60° , i = 75° and i = 90° . The minimum 
inclination angle of i = 45° is calculated using the relation 
between apparent flattening, g, intrinsic flattening, ^, and 



3 CONSTRAINTS ON THE DISTRIBUTION 
OF STELLAR VELOCITIES 

As kinematic constraints, we combine long-slit kinematics 
probing the central region and planetary nebula (PN) LOS 
velocities that probe the outer regions of NGC 4649. Figure 
[2] illustrates the spatial coverage of the kinematic data. We 
assume that the x-axis increases along East and that the 
y-axis increases along North, and then we rotate the coordi- 
nate system to align the a;-axis with the major axis and the 
j/-axis with the minor axis. 



^ i = 0° corresponds to a face-on system, and i = 90° corre- 
sponds to an edge-on system. 
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3.1 Long-slit kinematics 
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Figure 3. Comparison of PN data with stellar surface-brightness 
and long-sht data, along the major axis: ( a) The surface- 
bright ness profile measured from photometry in lKormendv et alj 
(|2009|) (black dia monds) and the scaled n umber density calcu- 
lated from PNe in iTeodorescu et alJ 1 I2OIII) (black circles), along 
with Poissonian errors, (b) and (c) show mean velocity and ve- 
lo city dispersion pro files respectively, measured from long slits 
in lPinknev et aU 2003) (black diamonds) , adapted from long-slit 
kinematic data in De Bruvne et al"] 1 I2OOII ) (red diamonds, not fit- 
ted in the NMAGIC models), and from PNe (black circles). 



We use the Ions -slit absorption-line kinematics of 
IPinknev et"al] (j2003l ). obtained from STIS/HST measure- 
ments and ground-based spectroscopic measurements on tlie 
2.4m MDM telescope. They were measured along position 
angles of 105° (the major axis), 127°, 133° and 173° mea- 
sured from North towards East in the sky, and extending to 
69", 28", 44", and 44" respectively. The lo cation and ori- 
entation of these slits are shown in Figure iPinknev" et al.l 
(2003) derive profiles of the line-of-sight (LOS) mean veloc- 
ity, velocity dispersion and the higher-order Gauss-Hermite 
moments and /i4. Figures [S^b) and (c) show the mean ve- 
locity and velocity dispersion profiles measured by the long- 
slit data along the major axis. The mean velocity increases 
from zero at the centre to about 110 km/s at ~ 45" and then 
decreases to about 95 km/s at ~ 70". The velocity disper- 
sion decreases from about 410 km/s at the centre to about 
270 km/s at ~ 70". For comparison we have overplotted 
major-axis long-slit kinematic data from iDe Bruvne et all 
(2001), averaged over both sides o f the galaxy. The agree- 
ment is very good. The data from IPe Bruvne et al.l ()200ll ) 
are not used in the NMAGIC modelling. 



3.2 Planetary nebulae 

To enable us to probe the mass and orbital structure in the 
halo of NGC 4649, we use LOS velocities derived from obser- 
vations of PNe with FORS2 on the VLT and FOCAS on Sub- 
aru. The observations and calculat ions of the velocities ar e 
described in a companion paper (jTeodorescu et al] I2OIII ) . 
Contaminants from the neighbouring spiral, NGC 4647, are 
remov ed using the technique developed by iMcNeil et al.l 
(|2010l ). This results in a catalogue of 298 PNe extending 
from 27" to 416", therefore overlapping with the long-slit 
kinematic data. 

The location of the PNe are shown in Figure [2] and Fig- 
ure [3J a) compares the measured major-axis V^-band surface 
brightness profile with the scaled number density profile of 
the PNe. The PN number density profile is calculated in a 
cone of angular width 30° centred on the major axis, and 
then scaled to match the photometry. We do not consider in- 
completeness corrections, which are especially important in 
the centra l region, where the br ight stellar component masks 
PNe (e.g. ICoccato et aljlioogi ). However, the PNe number 
density and surface-brightness profiles agree well, even at 
the innermost point. 

Figure [3l^b) shows the mean velocity profile of the PNe 
calculated in a 30°-cone centred on the major axis. In the 
region of overlap with the long-slit kinematics, the mean ve- 
locity profiles agree and further out the PN mean velocity 
profile decreases outwards to about -20 km/s at ~ 310". 
The large errors bars associated with the PN mean veloc- 
ity profile allow for a much shallower decrease in the outer 
parts, consistent with zero at 310". The combined long-slit 
and PN major-axis profile shows that the rotation of the 
stars increases to about 100 km/s at 50" and then starts 
decreasing outwards. 

Figure [Sjc) shows the velocity dispersion profile of the 
PNe calculated in a 30°-cone centred on the major axis. 
The profile agrees with the long-slit kinematic profiles in the 
region of overlap and then continues to fall less steeply to a 
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Figure 5. Minor-axis mean velocity profile: Tlie filled, black cir- 
cles show the mean velocity profile derived from PNe along a 
30°-cone centred on the minor axis. The open, black circles show 
the mean velocity profile derived from PNe along a pseudo slit 
of width 60" centred on the minor axis. The solid black bars 
denote the l-tr errors and the dotted black bars denote the 2- 
a errors, estimated from generating pseudo data sets. The solid 
black line shows the mean velocity profile extracted along the 
minor axis of the smoothed mean velocity field derived from the 
PNe. Red diamonds sh ow long-slit kinematic data adapted from 
iDe Bruvne et al.l |20o3). 



value of about 200 km/s at ~ 310". The consistent stellar 
and PN density and kinematic profiles imply that PNe are 
good tracers of the stars both in density and kinematics. 
Figures |4j a) and (b) show smoothed mean velocity and 



velocity dispersion maps for the region of NGC 4649 covered 
by t he PNe, created by th e PNpack suite of IDL routines 
from ICoccato et al] (|2009l ). Rotation approximately about 
the minor axis is seen clearly in (a). There is also a group 
of PNe in the top of the field that suggest some rotation 
about the major axis, a signature of triaxiality. To inves- 
tigate the strength of this signature, we first extract mean 
velocity profiles for the PNe along a pseudo slit of width 60" 
and a cone of angular width 30°, centred on the minor axis. 
Then to estimate the errors on these profiles, we generate 
100 pseudo sets of PN LOS velocities at the same positions 
as the original data, assuming that the line-of-sight velocity 
distribution (LOSVD) at each point in space is given by the 
mean velocity and velocity dispersion of the smoothed kine- 
matic fields. For each of the pseudo data sets, we extract new 
mean velocity profiles as done for the original data set. l-cr 
and 2-G errors are estimated from the spread in the mean 
velocity values at each position. Figure [5] shows the result- 
ing minor-axis mean velocity profiles along with the errors. 
We also extract a mean velocity profile along the minor axis 
from the smoothed mean velocity field. The profile along the 
pseudo slit is almost consistent with zero within 1-a errors 
and both the pseudo-slit and cone profiles are consistent 
with zero within 2-a errors. The profile derived from the 
smoothed mean-velocity map is very similar within ~ 90" 
but then decreases to zero, therefore suggesting a lower mag- 
nitude of velocity compared to the other two profiles. This 
profile should however be treated with caution in the outer 
parts as it is based on far fewer PN e there. The long-slit kine- 
matics of iDe Bruvne et al.l (|200ll ) also show some rotation 
along the minor axis at about 20 km/s between 50-80". We 
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do not plot their last point of ~ 90 km/s at the last radius 
of ^ 95", as it appears unphysically higher than the rest of 
the profile, ft would however still be in agreement with the 
PN minor-axis mean velocity profiles and is much less than 
the measured velocity dispersion. Therefore there is some 
evidence for triaxiality but it is weak, making it difficult to 
constrain the viewing angles of the system. Therefore we 
will proceed with models assuming an oblate, axisymmetric 
stellar distribution. 



4 DYNAMICAL MODELS 

Here we describe how we set up the initial model and how we 
prepare the photometric and kinematic target observables 
for creating dynamical models with NMAGIC. We then de- 
scribe how we obtain models first fitting to the light and 
long-slit constraints, and then models fitting light, long-slit 
and planetary nebula (PN) constraints. 

4.1 NMAGIC 

NMAGfC is an N-body made-to-measure code described in 
Ide Lorenzi et al] l|2007l ) that has been successfully applied 
to the intermediate-lum inosity elliptical ga laxies NGC 4697 
Jde Lorenzi et al.l 12008^) and NGC 3379 (|de Lorenzi et all 
l2009l ). ft finds the best intrinsic distribution of stars and 
their velocities that projects to fit the photometric and kine- 
matic data. The co de builds on the particle- based made-to- 
measure method of ISver fc Tremaind l| 19961 ) by accounting 
for observational errors and therefore allowing for an as- 
sessment of the quality of the model fits to the target data. 
T here have also b een r ecent implementatio ns of this method 
bv lDehnenI (|2009l ) and lLong fc Mad (|2010l ). 

NMAGIC starts with some initial particle model, where 
each particle has 3-D spatial coordinates, 3-D velocity co- 
ordinates and a weight. The particles are advanced accord- 
ing to the gravitational force to sample their orbits and the 
weights of the particles are adjusted simultaneously accord- 
ing to the force-of-cha nge (FOG) eq uation, given in Equa- 
tion (10) of [d o Lorenz i etal] l|2007l 'l and Equation (22) of 
Ide Lorenzi et al.. (,2008 ). This equation maximises the merit 
function, F defined as: 

F^^JiS-\x^+C (2) 

where S is a profit function, measures the goodness of 
fit to the density and long-slit target observables, and C is 
a likelihood term measuring the goodness of fit to the PN 
target observables. The parameter ^ changes the balance 
between the profit function and the goodness-of-fit terms. 
For the profit function, the entropy of the weights is used: 

^--E-i<5) (3) 

where Wi are the initial or prior weights of the particles, and 
Wi axe the current weights of the particles. The entropy term 
causes particle weights to remain close to their priors. There- 
fore a higher value of will generally result in smoother 
intrinsic distributions. We have chosen a value of 2 x 10"^, 
based on p revious work using NMAGIC (|de Lorenzi et all 
I2008l . l2009h . 
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Figure 6. Circular velocity curves of NGC 4649: The solid 
black line ( VCl) shows the best-fit circular velocity curve from 
iDas et al.1 l l2010r) . The solid re d line (VC2) sho w s the best- 
fit circular velocity curve from IShen fc Gebhard 3 ll20ld ). The 
blue lines show the circu lar velocity curves based on that from 
IShen fc Gebhardtl l|2010h but with the outer slope reduced to 
0.106 (VC3, dotted), 0.000 [VCi, solid) and -0.106 iVC5, 
dashed) . The green line show s the circular velocity curve based 
on that from lDas et al.l 1I2OIOI ) but with the outer slope reduced 
to 0.000 {VC6). The vertical black dash-dotted line shows the 
maximum radial extent of the kinematic constrain ts. For com- 
parison, the circular velocity curves o btained by Hu mphrey et all 
feooa) (soUd. pink).lHumphrev et al.l 12008 ) (dashed, pink), and 
iNagino fc Matsushital ||2009|) (cyan) are overplotted. 

4-. 1.1 The gravitational force 

In this work we treat the whole gravitational force as an 
external force acting on the system, and therefore it can 
take any form. In the first instance we co nsider the c i rcular 
velocity curves derived for NGC 4649 b y iDas et al.1 (|2010|) 
(VCl) using X-ray observations and IShen fc Gebhardtl 
(j201Q) {VC2) using optical observations and dynamical 
models. These two mass distributions only differ in the cen- 
tral ~ 12 kpc and are shown in Figure |6] The X-ray observa- 
tions give information on the temperature and density pro- 
files of the hot gas in massive elliptical galaxies. If the hot 
gas is approximately spherical and in hydrostatic equilib- 
rium, then the temperature and density profiles are related 
to the circular velocity curve by: 

^.^fc^dlnP 
finip d In r 

where T is the temperature of the gas and r is the 3-D radius 
from the centre of the galaxy, jj, = 0.61 is the average gas 
particle mass in terms of the proton mass, rrip. This value of 
fi corresponds to a helium number density of 7.92 x lO"'^ and 
0.5 solar abundance of heavier elements. We assume that the 
gas is ideal and therefore the gas pressure P = nksT, where 
n is the particle number density of the gas. iDas et al] (|2010l ) 
applied hydrostatic equilibrium using a new non-parametric 
Bayesian approach to density and temperature profiles of 
the hot gas d erived from Chan dra and XMM-Newton obser- 
vations (.Churazov et alj|2010l ). 

The best-fit mass profile of IShen fc Gebhardtl (|2010l ) 
was obtained from axisymmetric Schw arzschild models 
assuming a stellar density profile from iKormendv et al] 
l|2009l ) and fitting to long-slit kinematic constraints 
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from IPinknev et al.l l|2003l ') and GC LOS velocities from 
iHwang et alj l|2008l ). 

In both cases we assume that the gravitational force is 
spherically symmetric. Therefore as the stars are assumed to 
have an oblate, axisymmetric distribution, the dark matter 
halo must be prolate in the centre where the stars domi- 
nate, and almost spherical in the outer parts where the dark 
matter dominates. 



4.1.2 The initial particle model 

We set up an initial model of 750000 particles extending to 
~ 205 kpc. We assume a density distribution of the par- 
ticles given by a spherical deprojection of the circularly- 
averaged surface-brightness profile. We define the intrin- 
sic stell ar velocity distr ibution using the circularity func- 
tions of iGerhardI (|l99ll ). resulting in an anisotropy profile 
(/3 = 1 — ag/a^) that is isotropic in the centre but mod- 
erately radial (/? ~ 0.5) in the outer parts. T his choice re- 
flects the r esults of num e rical simulations by lAbadi et al.l 
l|2006t ) and lOfiorbe et al.l (|2007l ). We then solve for the en- 
ergy part of the distribution function assuming a poten- 
tial given by VCl. The particles' coordinates and velocities 
are c hosen according to the c omple te distribution function 
after iDebattista fc SellwoodI (|200d ) and they are assigned 
equal weights of Wi = 1/750000 to produce the initial par- 
ticle model. As the gravitational field is fixed with time in 
our models, the orbits of the particles do not change with 
time. Therefore our initial particle model is conceptually 
very similar to the orbit libraries used as initial conditions 
in Schwarzschild models. 



4.1.3 Preparation of stellar density target observables 

We represent the 3-D luminosity density profiles obtained 
for each inclination from the deprojection in terms of Aim 
coefficients (dc Lorcnzi ct al. 2007 ) between 0.3 kpc and 201 
kpc. 0.3 kpc is the innermost radius of the X-ray circular ve- 
locity curve (VCl) and 201 kpc is the radius at which the 
density of the initial particle model starts falling off from 
the deprojected density profile. As the stellar distribution 
is assumed to be axisymmetric we set all moments that de- 
scribe non-axisymmetry to zero. We calculate the Ai^s over 
a grid of 60 radii. Poissonian errors are assumed for the ra- 
dial mass profile and Monte-Carlo simulations are used to 
determine errors for the higher-order mass moments. 



4.1.4 Preparation of kinematic target observables 

As we are assuming that the stellar distribution in NGC 
4649 is an oblate, axisymmetric system rotating about its 
minor axis, for every mass element at {x,y,v), we can as- 
sume that there is also a mass element at {~x, ~y, —v) and 
a mass element at {—x,y, ~v). Therefore we can add three 
more long slits at position angles of 105° — 68° — 37°, 
105° - 28° = 77° and 105° - 22° = 83°, which are spa- 
tial reflections around the x-axis of the long slits positioned 
at 173°, 133° and 127°. The PNe are increased 4-fold by 
imposing the above reflections, resulting in a sample of 1192 
PNe. As NMAGIC is a particle-based code, the weight of 



in 


^alms 


-^long-slit 




-F 


(1) 


(2) 


(3) 


(4) 


(5) 


45 
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0.632,0.521 
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644.6,641.6 
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0.220,0.196 


0.655,0.518 


0.396,0.326 


761.5,661.5 



Table 1. Best-fit NMAGIC models for NGC 4649 assuming dif- 
ferent inclinations in the X-ray {VCl, left in columns (2)— (5)) 
and dynamical {VC2, right in columns (2)-(5)) potentials : (1) 
Inclination assumed for model, (2) per target observable, 
(3) psr long-slit target observable, (4) P^'' target observ- 
able, and (5) the merit function. There are 960 density target 
observables and 655 long-slit target observables. 




Figure 7. Fit to the first moment of the Ai^s or differential 
stellar mass distribution for an inclination of 75° in the circular 
velocity curves VCl (black) and VC2 (red). 



each particle will be changed more evenly over each orbit 
with the 4-fold sample compared to the unfolded sample. 

For NMAGIC the light in each of the slit cells needs to 
be calculated, as the light- weighted kinematics are fit rather 
than the kinematics directly. The light in the cells is calcu- 
lated by integrating the surface-brightness distribution over 
the dimensions of each of the slit cells (assumed to have a 
width of 5") using a Monte-Carlo integration scheme. 

To fit the P N LOS velocities , we u s e the likelihood 
method used in Ide Lorenzi et al. I l|2008l . I2OO9I ). In this 
method the particles are binned in radial and angular seg- 
ments to calculate the LOSVD in each segment. Then the 
likelihood of each PN belonging to the LOSVD of its segment 
is maximised for all the PNe. For the binning, we choose 3 
radial bins (corrected for an average projected ellipticity of 
0.857) and 6 equally-spaced angular bins with the first cen- 
tred on the major axis. 



4.1.5 Logistics of a 'run' 

For each run we fit the desired data starting with the initial 
particle model generated above. There is an initial relaxation 
phase of 2000 steps where the particles are advanced accord- 
ing to the force and their weights are not changed. Each step 
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Figure 8. Fits to the long-slit data (a) along the major axis and (b) along a slit placed at a position angle of 133°: From the top down 
are shown mean velocity, velocity dispersion and the Gauss-Hermite moments h-j, and /i4. Black circles show the data, the black region 
shows the models for the various inclinations probed in VCl, and the red region shows the models probed for the various inclinations in 
VC2. 



is equivalent to 4.4 x 10^ yrslf] We then have another 2000 
steps during which we temporally smooth the observables 
and the weights of the particles are still not changed. Then 
the core phase of the run starts, where the particles are ad- 
vanced according to the gravitational force and the weights 
of the particles are changed according to the FOC equation. 
Once the fractional change in the total fa-Us to less than 
1.5xlO~^ over 5000 steps, we define the model to have con- 
verged. Finally this is followed by a free evolution, where 
the particles are advanced for a further 20000 steps and the 
weights are no longer changed, to check that the converged 
model has been sufficiently phase-mixed. The particles are 
advanced using an adaptive leapfrog scheme. 

In the FOC equation, the respective contributions of 
the Aim, long-slit kinematic and PN kinematic terms can 
be very different. The errors are very small for the AimS 
but much larger for the PNe, resulting in a much larger 
contribution from the Aim terms to the FOC equation. As 
a result, changes made to the weights of the particles by the 
PN data will be small, and need to be made many times 
before a converged model is reached. Therefore as the PNe 
and to some extent the long-slit are important in the halo, 
after some tests, we increase their contributions to the FOC 
equation by factors of 20 and 2 respectively, to ensure that 
the halo is adequately modelled. 



mass model of IShen fc Gebhardtl (|2010l ') {VC2) and the X- 
ray mass profiles in the literature, which are lower in the 
central ~ 12 kpc. F or the X-ray mass profile, we choose that 



oflDas et all ||2010|) (VC l), which is very similar to that of 
IShen fc Gebhardtl (|2010l ') in the outer parts. As the poten- 
tials only differ in the central 12 kpc, and the PN kinematic 
constraints are much weaker than the density and long-slit 
kinematic constraints, we will not include the PNe for these 
models. This saves considerable computational time because 
the PNe are further out and therefore the particles need to 
be integrated for longer there to achieve convergence. Addi- 
tionally the likelihood method is itself computationally ex- 
pensive. 

We carry out NMAGIC models assuming inclinations of 
45°, 60°, 75° and 90° for the stellar distribution, for which 
we also carried out the deprojections of the photometry. Ta- 
ble [T] shows the P^r data point for the density observ- 
ables, the long-slit kinematics, all the observables, and the 
merit function F of the final models in each inclination and 
potential. NMAGIC fits to the light-weighted kinematic ob- 
servables as well as the light in each cell, and therefore the 
long-slit is calculated as a sum over these variables. 

We find that both potentials prefer (in the sense of a 
lower x^) a-n inclination of 75° of the stellar system. In this 
inclination, neither potential is preferred by the combination 
of the photometric and long-slit constraints. 



4.2 Models fitting density and long-slit kinematic 
constraints only 

First we would like to understand whether the observational 
constraints are able to differentiate between the dynamical 



^ As we are approaching equilibrium systems, the physical time 
that lapses is less important than the number of steps required 
to achieve it. 



4-2.1 Fits to the observables 

Figure [7] shows the fits in both potentials to the first mo- 
ment of the AimS or differential stellar mass distribution 
for the most favoured inclination of 75°, in potentials VCl 
and VC2. The fits are almost indistinguishable from each 
other, and match the density constraints very well. Figure 
[5] shows the fits to the v, a, /13 and /14 moments along the 
major-axis slit and the slit placed along 133° , for both poten- 
tials and all inclinations. In general the fits are very similar 
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Potential 


Slit PA (°) 


V 


a 


ha 


hi 


X-ray 


105 


0.131 


-0.196 


-0.013 


0.032 


(VCl) 


127 


-0.005 


-0.048 


-0.036 


-0.179 




133 


0.110 


-0.015 


-0.014 


0.165 




173 


0.041 


-0.170 


0.033 


-0.003 




All 


0.073 


-0.111 


-0.007 


0.009 


Dynamical 


105 


0.041 


0.113 


-0.012 


0.149 


iVC2) 


127 


-0.043 


0.086 


-0.013 


-0.141 




133 


-0.032 


0.228 


0.018 


0.188 




173 


0.038 


0.128 


0.031 


0.067 




All 


0.019 


0.139 


0.006 


0.073 



2.5 



Table 2. Systematic differences between model and observed v, 
a, hs and h4, averaged over each of the first four slits for radii 
outside 4" , and then averaged over all the slits. 



but there are some systematic differences. The models in 
VCl produce mean velocity profiles with a lower magnitude 
(~ 5-20 km/s), a lower velocity dispersion (~ 10-30 km/s), 
an /i3 moment with a lower magnitude (0.01-0.02) and an 
/i4 moment that is very marginally lower on average. These 
differences are generally a small fraction of the error bars. 

All the models find a per data point less than 1. 
Then if the number of degrees of freedom is approximately 
equal to the number of data points, then we can be satisfied 
that we are fitting the data well, and even over-fitting. In 
reality however, the number of degrees of freedom is difficult 
to estimate because it is equal to the number of constraints 
subtracted by the number of parameters. The number of 
parameters is equal to the number of model parameters (e.g. 
halo, M/L, inclination) plus the number of weights that we 
are fitting, which is the number of particles. The number 
of constraints is equal to the number of data points plus 
the number of constraints introduced by the profit function, 
which is difficult to quantify. 

We now attempt a Ax^ analysis as done for example in 



IShen fc Gebhardtl(|20ld ). If we say each of our models could 
be characterised by four parameters (inclination, mass-to- 
light ratio and two parameters for the dark matter halo) that 
need to be d etermined, then A y^ for a 68.3% confidence in- 
terval is 4.7 (|Press et al.lll98^ ). or 0.003 per data point for 
1615 observables (960 Ai^ and 655 long-slit target observ- 
ables). Therefore in the absence of systematic errors, any 
models with a of more than 0.003 per data point greater 
than the minimum per data point can be ruled out be- 
cause we know they must definitely be outside the 68.3% 
confidence range around the true best model. Table [T] shows 
that the density observables most prefer an inclination of 
75° but the long-slit observables most prefer an inclination 
of 60°, in both potentials. Considering the data altogether, 
both potentials prefer an inclination of 75°. Generally the 
dynamical potential is preferred over the X-ray potential ex- 
cept for in the most favoured inclination, where VCl and 
VC2 are equally preferred. As the remaining models achieve 
a of more than 0.003 per data point away from the min- 
imum x^ of 0.304 per data point, we would rule them out 
with 68.3% confidence using the Ax^ approach. 

We also quantify the systematic difference, AS, between 
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Figure 9. Intrinsic kinematics: The left set of panels show the 
ratio of the radial velocity dispersion to the tangential velocity 
dispersion (bottom), and the ratio of the azimuthal velocity dis- 
persion to the meridional velocity dispersion (top), along R. The 
right set of panels show the same along z. The black region corre- 
sponds to the models carried out for various inclinations in VCl 
and the red region corresponds to the models carried out for var- 
ious inclinations in VC2, fitting photometric and long-slit con- 
straints only. The dashed black line shows the radial extent of 
the kinematic constraints. 



the model and data in both potentials: 



As = -y 



K 



M,j 



Ko,j 



(5) 



where K represents the observable v, a, or hi, M is the 
model value, O is the observed value, n is the number of 
observables, and e^.j is the error on the observed value. 
AS is calculated for each moment averaged over each of the 
first four slits (the remaining three slits are only refiections 
of the first three) for radii outside 4", and then averaged 
over all the slits. Only the results in the most favoured in- 
clination of 75° are shown in Table (5] In VCl, the mag- 
nitude of the model mean velocity is a little higher, while 
the velocity dispersion is systematically a little lower than 
the observations. The systematic differences in the /13 and 
hi moments are smaller. In VC2, the systematic differences 
are biggest also in the velocity dispersion, which is higher 
than the observations on average, and in the /14, which are 
also systematically higher. Systematic differences between 
the model and the data seem similar in both potentials and 
are also comparable to deviations used to rule out models 
(y/Ax^ = 0.003 = 0.055 per data point). This implies that 
neither mass model is exactly correct and that the true mass 
distribution lies somewhere in between, or that systematic 
effects (e.g. triaxiality) play a role. Therefore the Ax^ ap- 
proach should be used with caution. 



4-2.2 The orbital structure 

As we have assumed an oblate, axisymmetric system, we can 
average the intrinsic kinematics over the azimuthal angle, 
<j). This reduces the spatial coordinates to R in the equa- 
torial plane, which is the major axis of the system, and z 
in the meridional plane, the minor axis of the system. Fig- 
ure [9] shows the ratio of radial to tangential velocity dis- 
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VC2 


0.268 


0.510 


0.366 


3107.0 


4232.3 


VC3 


0.201 


0.496 


0.321 


3097.1 


4087.1 


VCi 


0.196 


0.489 


0.315 


3087.0 


3988.7 


VCh 


0.236 


0.480 


0.335 


3079.7 


4006.0 


VC6 


0.213 


0.818 


0.458 


3067.2 


4091.3 



Table 3. Best-fit NMAGIC models for NGC 4649 for an inclina- 
tion of 75° in a range of circular velocity curves: (1) Circular ve- 
locity curve, (2) per target observable, (3) per long-slit 
target observable, (4) P^r target observable, (5) log-likelihood 
of PNe belonging to particle LOSVDs, and (6) the merit function. 



0.3 - 




t 

50 100 150 200 

r (kpc) 

Figure 10. Fit to the first moment of the Aj^s or differential 
stellar mass distribution for an inclination of 75° in potentials 
VC2 (solid, red), VC3 (dotted, blue), VCi (solid, blue), VC5 
(dashed, blue), and VCd (green). 

persions (ar/crt), and the ratio of azimuthal to meridional 
velocity dispersions {o^/ag) along R and z. The tangen- 
tial velocity dispersion is defined as at = y^(oJ~fo^)/2. 
Within the radial extent of the kinematic constraints, there 
is a bias towards radial orbits in VCl along R and z with 
{arlot)max ~ 1-5. In the same region the orbital structure 
in VC2 is almost isotropic along R and z. On average the ra- 
tio of the radial to tangential velocity dispersions is 1.2-1.3 
higher in VCX than in VC2 throughout. The two compo- 
nents of tangential velocity dispersions have similar contri- 
butions in both potentials along z, and in the equatorial 
plane the azimuthal dispersions dominate. In both poten- 
tials, variation with inclination is greatest in the ratio be- 
tween the azimuthal and meridional dispersions along R. 

4.3 Models fitting density, long-slit kinematic 

constraints and planetary nebula line-of-sight 
velocities 

Now we will carry out models including also the PN kine- 
matics, so that we can check whether the PN kinematics 
are consis t ent w ith t he halo in the m ass determinations of 
iDas et all (|2010l ') and lShen fc Gebhard t (2010). As both po- 
tentials have almo st the same halo we first c arry out a model 
in the potential of lShen fc Gebhardd (|2010l ) (FC2). We find 
that the model velocity dispersions are systematically too 
high and therefore repeat models in a range of potentials 
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Figure 11. Fits to the long-slit data: From the top down are 
shown mean velocity, velocity dispersion and the Gauss-Hermite 
moments /13 and /i4 for an inclination of 75° and for the slit placed 
along the major axis. Black circles show the data and the lines 
show fits in potentials VC2 (soUd, red), VCi (dotted, blue), yC4 
(solid, blue), VC5 (dashed, blue), and VC6 (green). 

with less massive haloes (see Figure[6]and Table|3}. Incorpo- 
rating the PNe, the models prefer a halo with a significantly 
lower circular velocity of ~ 463 km/s compared t o a va lue 
of 570 km/s at 45 kpc found in lShen fc Gebhardt! (|20ld ). 

^.3.1 Fits to the observables 

Figures 1101 [TT1 and [12] show the fits to the first moment 
of the AimS, the fits to the Gauss-Hermite moments along 
the major axis and the fits to the PN LOS velocities, in 
VC2, for an inclination of 75°. Table [3] shows the statistics 
of the fits obtained for the new model incorporating the PN 
data. Even though visually the fits look very much the same 
in the Ai^s and long-slit, the values for the fit to the 
AimS and long-slit are slightly higher than before, though 
still much less than 1. This shows that in this potential, 
by trying to fit the PNe, the fits to the AimS and long-slit 
are slightly compromised. Looking at the PN kinematics, it 
appears that the mean velocity profiles are fit well but the 
velocity dispersions of the model are systematically higher 
than that measured by the PNe. This implies that the PNe 
are not consistent with the outer part of VC2, and therefore 
also VCl (both potentials agree very well outside ~ 12 kpc) . 

4-3.2 Fits in less massive haloes 

As the PN dispersions are systematically lower than pre- 
dicted by the models, we also investigate additional circular 
velocity curves with less massive haloes (Figure |S} . We fit 
a straight line to VC2 outside 7.6 kpc and obtain a slope 
of 0.212. We create three additional curves where the outer 
slopes are 0.106 {VC3), 0.0 (VCi) and -0.106 {VC5). Fi- 
nally we consider a curve that follows VCl until 4.9 kpc 
and then also has an outer slope of 0.0 {VC6). The model 
fits in these potentials are overplotted on Figures \TU\ [TT] [T^] 
and the statistics of the fits are given in TableO The x^ val- 
ues show that the AimS and long-slit data most prefer VC4, 
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Figure 12. Fits to tlie PN data for an inclination of 75°. Eacli of the plots in the bottom three rows shows the LOSVD in segments. 
Going upwards are segments at radii of 64", 145" and 316" and going towards the right are angular segments centred on 0° (major 
axis), 60°, 120°, 180°, 240°, and 300°. The kinematics along the latter three segments are reflections of the kinematics along the first 
three segments because of the oblate, axisymmetry imposed on the observational constraints. The top two sets of panels show the mean 
velocity and velocity dispersion profiles along the angular segments. The fits are in potentials VC2 (solid, red), VC3 (dotted, blue), VC4 
(solid, blue), VC5 (dashed, blue), and VC6 (green). 



and the likelihood values show that the PNe most prefer 
VC6. The merit function _F is a combination of the val- 
ues, likelihood and entropy, and this is a maximum in the 
potential 1/(74. Therefore we consider this the best of the 
potentials tried. The differences between the models in the 
AimS and fits to the long-slit kinematics are small except in 
a and /i4, which are systematically lower in VC6 compared 
to potentials VC2~VC5. Figure [T2l shows that the PN ve- 
locity dispersions are however very sensitive to the mass in 
the halo. 



4-3.3 The orbital structure 

Figure [13] shows that within ~ 6 kpc the intrinsic velocity 
dispersions are isotropic to mildly radial along R and z in 
potentials VC2-VC5 and moderately radial in VC6, con- 
sistent with the findings in Section [4.2.21 Further out until 
the radial extent of the PN kinematic constraints, the or- 
bital structure in all potentials becomes moderately radial 
along R and more isotropic along z. The two components of 
tangential velocity dispersions have similar contributions in 
all the potentials, and in the equatorial plane the azimuthal 
dispersions dominate over meridional velocity dispersions. 
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Figure 13. Intrinsic kinematics: The left set of panels show the 
ratio of the radial velocity dispersion to the tangential velocity 
dispersion (bottom), and the ratio of the azimuthal velocity dis- 
persion to the meridional velocity dispersion (top), along R. The 
right set of panels sho w the same along z. The d ash-dotted black 
line shows the results of lShen fc Gebhardtl 1I2OIOI ') who fit GC kine- 
matics instead of PN kinematics. The fits are done in potentials 
VC2 (solid, red), VC3 (dotted, blue), VCi (solid, blue), VC5 
(dashed, blue), and VC6 (green). The dashed black line shows 
the radial extent of the kinematic constraints. 
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5 DISCUSSION 

Here we discuss what we have learned about the dark mat- 
ter halo, orbital structure and inclination of the stellar sys- 
tem in NGC 4649. We also discuss whether the PNe and 
GC systems are dynamically consistent with each other, and 
whether the X-ray mass distributions can be used to deter- 
mine dark matter mass fractions and orbital structures in 
massive elliptical galaxies. 



5.1 The dark matter halo of NGC 4649 

Figure [6] shows the circular velocity curves that we have 
probed in this work, compared to other recent determi- 
nations from Chandra observations l|Humphrev et al.ll2006l . 
I2OO8I ) and a combination of C handra and XMM-Newton ob- 
servations (|Nagino fc Matsushita, 2009. '). 

Our models fitting only density and long-slit kinematic 
constraint s cannot d i stingu ish be tween the circular v elocit y 
curves of lOas et all (|2010[ ) and IShen fc Gebhardtl (|2010h . 
which differ only in the central ~ 12 kpc. A look at the 
systematic differences between the models and observations 
in these potentials suggests that the true circular velocity 
curve in the central ~ 12 kpc probably lies somewhere be- 
tween 425-500 km/s. 

Mo dels created incorporat ing additionally PN kine- 
matics (|Teodorescu et aTl I2OIII I prefer a circular velocity 
curve that is flat outside ~ 12 kpc at a value of ~ 463 
km/s. This is mos t consistent with the X-ray determina- 
tion of iNagino fc^M atsu shit a ( 2009 ), slightly higher than 
the determinations of iHumphrev et all (|2006| . l200a) . and 

and 
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lower than the determ inations of 
IShen fc Gebhardtl (|2010t ). The sensitivity of the PN velocity 
dispersions to the circular velocity curve in the halo may be 
a result of a light density profile that falls off appr oximately 
as -3 and an almost flat circular velocity curve. I Gerhard 
l|l993l ) showed that for such systems, the projected velocity 
dispersions are constant and independent of anisotropy. The 
constant circular velocity Vc is then related to the constant 
projected velocity dispersion ap by Vc — x ap. The av- 
erage velocity dispersion of the PNe outside ~ 70" is ~ 240 
km/s, therefore predicting a circular velocity in the halo of 
416 km/s. This lies in between the value of 383 km/s of VC6, 
the halo most preferred by just the PNe, and the value of 
463 km/s of VCA, the halo most preferred considering all 
the data together. 

To estimate the dark matter mass fraction, we first 
calculate the luminosity enclosed within radius r by inte- 
grating over the luminosity density profile obtained from 
the spherical deprojection of the surface-brightness pro- 
file. To obtain the mass in stars we multiply this with a 
cons tant stellar population mass-to-light ratio estimated in 
iTeodorescu et a l. (2011). They used ages and metallicities 
in lTrager et al.l ( pOOCj) and the evolutio nary population syn- 
thesis models of lMarastonI (|l998l . [200^ for a Salpeter initial 
mass function with a lower mass limit of O.IMq. Their value 
of 10 in the B band corresponds to a value of about 7.6 in the 
V band. Assuming that the mass in gas is negligible, we sub- 
tract the mass in stars from the total mass. Figure [141 shows 
the dark matter mass fractions corresponding to the poten- 
tials we have explored and the potentials from the literature 
shown in Figure [G] The dark matter mass fraction corre- 
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Figure 14. Dark matter mass fraction of NGC 4649 ac- 
cording to potentials VCl (black), VC2 (red), VC3 (solid, 
blue), VCi (dotted, blue), VC5 (dashed, blue), and VC6 
(green). For comparison, the dark matter mas s fractions calcu- 
lated from circular velocity curves obtained by [Humphrey et all 
ll2006l'l (solid. pinkLlHumplirev et all \200^ (dashed, pink), and 
iNagino fc Matsushital 1I2OO9I ) (cyan) are also overplottcd. 



spending to our best potential VCA is ~ 0.39 at 10.5 kpc 
(IRe), 0.5 at about 16 kpc (1.5i?e) and ~ 0.78 at the radius 
of the last PN (~45 kpc ^ ARe). Figur elTilshows that our de- 
termination is most similar to that of lNagino fc Matsushital 
(200^ as for the circular velocity curve. 

There are also several determinations of dark mat- 
ter fractions for othe r ellip tical galaxies in the literature. 
INagino fc Matsushital l|2009l ) analysed Chandra and XMM- 
Newton observations for a sample of 22 elliptical galaxies 
and found equality between dark matter and luminous mat- 
ter at around 3Re and a dark ma t ter ma ss fra ction of around 
0.66 a t ~ 6Re. ICerhard et al] l|200lh and iThomas et ai] 
(|2007l ) found dark matter mass fractions within IRe of 10- 
40% and 10-50% respectively, for samples of nearb y and 
Coma cluster elliptical galaxies. ICerhard et al.l l|200ll ) found 
equality between the mass contributions of the dark mat- 
ter and luminous components at 2-4_Re, and for the Coma 
galaxies with sufficiently spati ally extended d a ta (ge nerally 
the less luminous ellipticals), iThomas" et all (|2007l ) found 
dar k matter mass fractions be tween 65-75% at 4i?e. 

iTreu fc Koopmani (|2004l ) used a combined lensing and 
stellar dynamical approach to obtain the dark matter mass 
fractions of elliptical and lenticular galaxies up to a redshift 
of 1, and found projected values of 0.15-0.65 within a cylin- 
der of radius IRe - For massive early-typ e galaxies taken from 
the Sloan Lens ACS S urvey (SLAGS), lAuger et all (|2009l ) 
and I Auger et all (|2010l ) found that the dark matter fraction 
within Re/2 ranges between ~ 0.3-0.7 assuming a Chabrier 
stellar IMF (the Salpeter stellar IMF gives lower dark mat- 
ter mass fractions in their study that are sometimes nega- 
tive), with more massi v e gala xies having higher dark matter 
mass fractions. iGrilld 1I2OIOI ) used simple mass models to 
estimate a dark matter mass fraction of 0.64o if projected 
inside a cylinder of radius IRe for a sample of approximately 
170000 massive, elliptical galaxies observed in SDSS, from 
which the SLAGS sample is obtained. We estimate the effect 
of projection on the dark matter mass fraction by assuming 
an NFW density profile for the dark matter profile, with a 
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virial radius of 300 kpc and a concentration of 10 (typical for 
elliptical galaxies) . The ratio of the dark matter mass within 
a cylinder of length two virial radii along the LOS and ra- 
dius Rs, to that within a sphere of radius IRe is about 2. 
Assuming that the ratio of the stellar mass between the two 
regions is almost 1, then the dark matter mass fraction cal- 
culated within a sphere of radius IRe would be lower, for 
exa mple 0.47 instead of 0.64. 

lOfiorbe et al.l (|2007l ') analysed the mass and velocity dis- 
tributions of elliptical-like objects at zero redshift in a set 
of self-consistent hydrodynamical simulations set in the cur- 
rent cosmological paradigm. They found that the objects are 
embedded in massive, dark matter haloes with dark matter 
mass fractions ranging between 0.3-0.6 at IRe- 

To summarise, there is a range in the dark mat- 
ter mass fractions at 0.5-1-Re in the literature, and the 
value we obtain for NGC 4649 is near the middle of this 
range. Fu rther out, the dark matter mass fraction s ob- 
tained bvlGerhard et"all ll200lh , iThomas et"all l|2007h and 
iNagino fc Matsushita (|2009l ') suggest on average a more dif- 
fuse dark matter halo than the one we find for NGC 4649. 
This may be because their samples include elliptical galaxies 
at a range of luminosities, while massive elliptical galaxies 
like NGC 4649 may have more massive dark m atter haloes 
(e.g. ICappellari et al.ll2006l : [Auger et allbOloD . 

5.2 Orbital structure in NGC 4649 

The central ~ 6 kpc of NGC 4649 has an isotropic {/3 — 
1 — (Jf/a^. ~ 0) to mildly radial (/3 ~ 0.4 ) orbital structure 
accord ing to the dynamical pot ential of Shen fc Gebhardtl 
l|2010l ) and the X-ray potential of iDas et al.l ( 201Cf ). between 
which we are unable to distinguish. If we assume that the 
true mass distribution in the central ~ 12 kpc lies some- 
where between these two, as suggested by the systematic 
differences between the models and the observations in the 
two potentials, then we can infer an orbital structure in this 
region of /3 ~ 0.2±0.2. Using additionally the PN con- 
straints and exploring several potentials, we find that the or- 
bital structure outside ~ 4 kpc becomes slightly more radial 
(/? ~ 0.5) along R, but more isotropic along z, with little de- 
pendence on the exact halo assumed. Along R, the azimuthal 
velocity dispersions are slightly higher than the meridional 
velocity dispersions throughout, indicating that the stel- 
lar system may be flattened b y a meridional anisotropy in 
the velocity dispers i on tensor ([Dehnen &: GerhardI Il993al lbl: 
iThomas et al. l l2009l ). lThomas" et al.l (|2009l ) also use axisym- 
metric toy models to show that flattening by meridional 
anisotropy maximises the entropy for a given density dis- 
tribution. Along z, the azimuthal and meridional velocity 
dispersions are equal, as one would expect for an oblate, 
axisymmetric system. 

There is no general consensus on the orbital structure 
in the outer parts of elliptical galaxies. Dynamical models 
fitting outer kinematics of intermediate-luminosity elliptical 
galaxies have found a mo derately radial orbital structure 
in th e halo of NGC 4697 (Ide Lorenzi et al.ll2008l ). In NGC 
3379, Ide Lorenzi et al] l|2009l ') found that systems with both 
isotropic and moderately ra dial orbital structures we re con- 
sistent with the data, and iNapolitano et all (|2009l ) found 
moderately radial anisotropy in the halo of NGC 4494. Dy- 
namical models fitting outer kinematics in more massive el- 



liptical galaxies such as NGC 1399 (|Schuberth et al.ll2010l . 
the GCs used are believed to trace the stellar kinema tics in 
this galaxy) and NGC 4374 l|Napolitano et al.1 12OI0I ') have 
found mildly radial and isotropic orbital structures in the 
halo respectively. For the elliptical galaxies in the Coma 
cluste r, which have a range of luminosities, IThomas et al.l 
(j2007h found mild radial anisotropy along the major axis 
and sometimes tangent ial anisotropy along the minor axis. 

From simulations, lAbadi et al.l ([2006) found that the 
outer haloes of massive ellipticals are strongly radial due to 
smaller galaxies tha t have been accreted on to the central ob- 
ject. The analysis of lOfiorbe et al.l l|2007l ) of elliptical-like ob- 
jects at zero redshift also found a radial orbital structure for 
the stars at an almos t constant value of /3 ~ 0.5 throughout. 
IThomas et al.l (|2009l ) analysed the orbital structure of colli - 
sionless disc merger remnants from iNaab fc BurkertI (I2003D 
and found them to be strongly radially anisotropic. 

To summarise, dynamical models show that the orbital 
structure in the halo of massive elliptical galaxies is isotropic 
to quite radial , but less so than expected from simulations. 
IThomas et al.l (|2009t ) arrived at a similar conclusion and sug- 
gested that this could be due to gas dissipational effects. 



5.3 The inclination of NGC 4649 

The surface-brightness and long-slit kinematic constraints in 
both potentials VCl and VC2 prefer an inclination of 75° 
out of the four inclinations we have probed. As we only ex- 
plore a coarse grid of inclinations at intervals of 15°, we can 
attach an error of 7.5° to our best inclination. Thus it ap- 
pears that even if the X-ray potential may not be completely 
correct throughout, it may be used to find an approximate 
value for the inclination of the system. Th is appears to be 
contra ry to the work done for example by iKrainovic et al.l 
(|2005h . who found a degeneracy in the determination of the 
inclination in their construction of axisymmetric models for 
the elliptical galaxy NGC 2974. However, we have assumed 
the same spherical total potential in all the inclinations for 
our first set of models, and only the stellar distribution was 
axisymmetric and varied according to the inclination. To 
understand this issue in more depth, a range of mass pro- 
files need to be explored in each inclination to see whether 
equally good but different mass profiles can be found for 
each of the inclinations. 



5.4 Are the PNe and GCs dynamically consistent 
with each other? 

The NMAGIC models show that the PN kinematics prefer 
a less massive halo and a less radial orbital structure than 
that derived from mo dels using t he stellar density profile 
and GC kinematics in lShen fc Geb hardt (2010). However if 
both the PN and GC systems are in equilibrium, the same 
halo should be preferred although the orbital structure may 
be different. To obtain a deeper insight into this discrepancy, 
we examine the photometric and kinematic constraints as- 
sumed for the PNe and GCs in more detail. Figure [I5ja) 
compares the velocity dispersion profile calculated in circu- 
lar rings from the PNe to t hose for red, blue, and all GCs 
calculated by iHwang et al. 

1 12008). The errors bars on the 
GC kinematics are not shown because they mask all the 
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Figure 15. Properties of tlie GC populations: (a) Velocity dispersion profiles of all GCs (filled, black diamonds), red GCs (filled, red 
diamonds) and blue GCs (filled, blue diamonds) compared to that measured by the PNe (open, black circles). The errors on the GC 
velocity dispersions are approximately twice as big as the errors shown on the PN velocity dispersions, (b) Surface number density 
profiles of all GCs (filled, black diamonds), red GCs (filled, red diamonds) and blue GCs (filled, blue diamonds) compared to the stellar 
surface- brightness distribution (solid, black line) and the surface number density of PNe (open, black circles). The GC points are not 
independent as they are obtained with a moving window average. 



other points but are approximately twice as big as the er- 
rors shown on the PN velocity dispersions. One can see that 
between 100" and 200" , the PN velocity dispersions appear 
to be in agreement with the dispersions measured by the 
blue GCs (and therefore the total sample because the blue 
GCs dominate in numbers). The velocity dispersion of the 
red GCs decreases to almost 150 km/s at ~210" and then 
increases dramatically to about 320 km/s at ~390". The last 
PN velocity dispersion point at ~350" calculated in a ring 
extending from 260-440" has a value of 203 ± 19 km/s. Av- 
eraging all the GCs velocity dispersions in this region gives 
~ 240 ±22 km /s for 61 GCs, wh i ch wo uld correspond to 
the values fit bv lShen fc Gebhardtl (|2010l 'l. This is only just 
consistent and therefore it is possible that the PNe and GCs 
trace different kinematics. 

Figure llSf b) compares the surface-brightness distribu- 
tion of the stars wi th scaled surface number densities calcu- 
lated from PNe in iTeodorescu et al.l (|201ll ) and GCs from 
iHwang et all l|2008r ). Between 200-450". the surface densi- 
ties of the whole GC population follows the stellar surface- 
brightness distribution well. Further in, the GC density pro- 
file is shallower than that of the stars and outside this re- 
gion the density profile falls off more steeply than the stars, 
though the associated errors are larger. Therefore the den- 
sity profile of the GCs appears to be different from that of 
the stars. 

The GCs could still be in equilibrium however but 
just form a separate dynamical system in the same poten- 
tial. This can be quantified using the spherical second-order 
Jeans equation, relating the second-order moments of the in- 
trinsic velocity distribution to the densi ty profile of the stars, 
and the potential in which they move |Binnev fc Tremainel 
1 19871 ): 

^{u{r)a^4r)) + ^^i.{r)aUr) + v(t)^ = (6) 
dr r r 

where Vc is the circular velocity curve, Or is the intrinsic 



velocity dispersion of the tracer in the radial direction, v is 
the number density of the tracer, and the anisotropy /3 = 

This equation shows that as the PNe and GCs are resid- 
ing in the same halo but have different density profiles and 
probably different kinematics, for the GCs to be in equi- 
librium, they must also have a different orbital structure. 
We can estimate the orbital structure most easily for the 
case of a power-law density profile, constant anisotropy and 
constant circular velocity curve. Therefore we fit a power 
law to the surface-brightness measured by the stars outside 
200" and find a best-fit index of nstars = —2.2. For the glob- 
ular clusters we find rigcs = — 3.0 (ignoring the second last 
point, which seems unphysically low). The power-law indices 
of the intrinsic density profiles are then -3.2 and -4.0 respec- 
tively. The projected velocity dispersion is ~ 240 km/s out- 
side 200", and we assume a circular velocity of 463 km/s 
in the halo from our best potential VC4. Solving the Jeans 
equation and the equation relating projected velocity dis- 
persions to intrinsic velocity dispersions, we find /9 ~ — 1 for 
the GCs. This is different from the average value o f /3 ~ 0.4 
outside 200" obtained bv lShen fc Gebhardtl (|2010D . who as- 
sume the stellar density profile for the GCs. Our value agrees 
better with the modest ly tangentiall y biased velocity ellip- 
soid inferred by Hwang et al.l l|2008l ) from spherical Jeans 
equations using the true G C density profile and t he mass 
distribution from X-rays in iHumphrev et al.l (|2006l ). which 
is more consistent with the mass distribution we find from 
the PNe. 

To summarise, we show that the GCs may be in equi- 
librium in the potential of NGC 4649, but this equilibrium 
is dynamically distinct from the stars and PNe. Therefore to 
infer the mass of the dark matter halo, the correct density 
profile and kinematics need to be used. 
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5.5 Are the mass distributions from X-rays 

accurate enough to determine dark matter 
mass fractions and orbital structures? 

Using the photometric and long- slit kinematic cons t raint s 
in the dynamical potential from Shen fc Gebhardtl (I2OI0I ) 
and the X-ray potential from lPas et al.l ( )201Cf ). we find that 
the same inclination of 75° is preferred for the stellar sys- 
tem. In this inclination we are unabl e to distinguish b e tween 
them . As the X-ray potentials fr om Humphrey et al.l l|2006l . 
l2008h and iNaeino fc Matsushital (|2009( l are similarly low in 
the central ~ 12 kpc, we would expect a similar result if 
one of these were used instead. Looking at the systematic 
differences between the models and observations in the two 
potentials suggests that the true mass distribution in the 
central ~ 12 kpc lies somewhere in between. By assuming 
the X-ray mass distribution, a maximum systematic error of 
0.4 in P is made, i.e. one finds a more radial orbital structure. 

From our results, we do not see compelling evidence 
for non-thermal pressure contributions in the gas or multi- 
phase components in the gas in the central ~ 12 kpc. If they 
exist, their effe cts are smaller than wou ld be inferred from 
the potential of lShen fc Gebhardd (l2010l) . This is consistent 
with the work of lBrighenti et al.l (|20o'9i rThev model the hot 
gas in NGC 4649 using 2-D gas-dynamical computations and 
find that a turbulent pressure is required in NGC 4649, but it 
has a much smaller contribution than the thermal pressure. 

In the halo, the mass distribution most preferred by in- 
corporating the PN d ata is less massiv e than that found in 
the X-ray analysis of IPas et al.l ll2010h. slightly mor e mas- 
sive c ompared to the halo found in ' Humphrev et al.l ||2006|. 
l2008f) . and similar to that found by Nagino fc Matsushital 
lj2009f ). The dark matter mass fractions derived in the halo 
depend highly on the mass distribution assumed. The or- 
bital structure in the halo however does not change much 
between the potentials explored. Still the important ques- 
tion arises: why do the X-ray mass distributions differ in 
the outer regions? Looking at the temperature and pressure 
profiles measured from Chandra and XMM-Newton obser- 
vations in Das et al.. (,2019' ). we find that at ~ 25 kpc, the 
temperature profiles are consistent but the pressure profiles 
are not. The pressure calculated from the Chandra data is 
higher than that calculated from the XMM-Newton data, re- 
sulting in a flatter pressure profile in the outer parts. Look- 
ing at Equation Q sh ows that this wil l result in a lower 
circular velocity curve. iDas et al.l l|2010l ) use both sets of 
observations but omit the final points due to uncertainties 
associated with the deprojec tion, and therefo r e do not use 
the Chandra point at 25 kpc.lHum phrev et al. I (2006) do use 
this point however. Therefore it seems that the outer slope 
of the mass distribution from X-rays can be quite uncertain 
and possible effects that need to be explored in more detail 
are deprojection issues, metallicity gradients in the hot gas, 
and outffows. 

To summarise, we believe that using the X-ray mass dis- 
tribution may lead to a systematically more radial orbital 
structure in the central region. Models with the X-ray mass 
profile however are able to derive the inclination of the stel- 
lar system and the orbital structure in the halo. Therefore 
until the uncertainties in the derivation of X-ray mass dis- 
tributions are better understood, it is best to use them in 
conjunction with a dynamical mass analysis. 



6 CONCLUSIONS 

We have created dynamical models of the Virgo elliptical 
galaxy NGC 4649, using the highly flexible made-to-measure 
N-body code, NMAGIC, and observational constraints given 
by surface-brightness data, long-slit kinematics, and plane- 
tary nebula (PN) velocities. We explore a range of poten- 
tials based on X-ray mass distributions in the literature, 
which are similar in the central regions, but have different 
outer slopes, and a dynamical potential derived from globu- 
lar cluster (GC) velocities and a stellar density profile. The 
GC dynamical model prefers more mass in the central region 
compared to the X-ray potentials, and is on the top end of 
the range of X-ray mass profiles further out. 

Our models are not able to differentiate between the X- 
ray and GC mass profiles in the central ~ 12 kpc, and the 
systematic differences suggest that the true circular velocity 
curve in the central ~ 12 kpc may lie somewhere between 
425-500 km/s. Therefore if non-thermal pressures or multi- 
temperature components exist in the central region, their 
contribution is less than previously inferred. 

Outside ~ 12 kpc, the observational constraints pre- 
fer a circular velocity curve that is fiat with a value of 
~ 463 km/s, m ost consistent with the X-ray determination 
of lNagino fc Mat s ushita (,2009. ) . The PN velocity dispersions 
are very sensitive to the circular velocity curve in the halo, 
possibly as a result of a stellar density profile that falls off ap- 
proximately as -3 and an almost flat circular velocity curve 
(Gerhard"l993') . The discrepancy between the halo mass pre- 
ferred by the PN kinematics and that corresponding to the 
GC dynamical model shows that if the GCs are in equilib- 
rium, they are dynamically distinct from the stars and PNe. 
Therefore the correct density profile and kinematics need to 
be used to infer the mass of the dark matter halo. 

We find a dark matter mass fraction of 0.39 at IRe 
for NGC 4649, which is generally in agreement with the 
values in the literature. At 4iie we obtain a dark matter 
mass fraction of 0.78, suggesting a more massive ha l o tha n 
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typical for the sampl es anal ysed by iGerhard et al] 
iThomas et al.l |20o3) and iNagino fc Matsushital 
This may be because their samples include elliptical galaxies 
at a range of luminosities, while massive elliptical galaxies 
like NGC 4649 may have more massive dark m atter haloes 
(e.g. ICappellari et al.ll2006l : fAuger et al.ll2010l ). 

We find an orbital structure that is isotropic to mildly 
radial in the central ~ 6 kpc, depending on the potential 
assumed. Further out, we find that the orbital structure be- 
comes slightly more radial along R, but more isotropic along 
z, with little dependence on the exact halo assumed. Along 
R, the azimuthal velocity dispersions are slightly higher than 
the meridional velocity dispersions throughout, indicating 
that the st ellar system may be flattened by a meridiona l 
anisotropy (jDehnen fc Gerhard|[l993al : iThomas et "aLll2009l ') . 
The orbital structure in the halo of NGC 4649 is less radial 
than expected f rom simulations, po ssibly due to gas dissi- 
pational effects (IThomas et al.|[2009l ). 

Assuming a mass distribution from X-rays leads to a 
systematically more radial orbital structure in the central 
region, but recovers the orbital structure in the halo. The 
inclination of the stellar system is also recovered. It is ap- 
parent that until the uncertainties in the derivation of X-ray 
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mass distributions are better understood, they should only 
be used alongside a more detailed dynamical mass analysis. 
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